Code
library (tidyverse)
library (magrittr)
library (gt)
library (patchwork)
library (SummarizedExperiment)
library (tidySummarizedExperiment)
library (brms)
library (gtsummary)
library (labelled)
library (scales)
# library(gridExtra)
library (ggpubr)
# library(mia) # BiocManager::install("mia")
theme_set (theme_light ())
tmp <- fs:: dir_map ("R" , source)
rm (tmp)
Code
tmp <- fs:: dir_map ("../VIBRANT-99-utils/R/" , source)
rm (tmp)
This document presents the primary analyses of the VIBRANT trial.
Loading the data
Code
mae_files <-
fs:: dir_ls (str_c (data_dir (), "04 unblinded MAEs/" ), regexp = "mae_full_.* \\ .rds$" ) |>
sort (decreasing = TRUE ) |>
magrittr:: extract (1 )
mae <- readRDS (mae_files)
rm (mae_files)
Analyses
Demographics (Table 1)
Summary table
We create a table1_data data.frame from the participant_crfs_merged table stored in the @metadata slot of the mae.
This table contains almost all variables needed to create table 1 (i.e., site, age, race, ethn, cut_size_meals, eat_less, hungry_did_not_eat, sexual_partners_lifetime, sexual_partners_past_month, contraception_method).
In addition, the mae@colData contains the participants’ arm and study population flags (here we keep the mITT population).
From the variables listed above, we create the variable food_insecurity using the three variables cut_size_meals, eat_less and hungry_did_not_eat: if any of these three variables is “Yes”, then the participant is considered to be affected by food insecurity.
Code
table1_data <-
metadata (mae)$ participant_crfs_merged |>
dplyr:: left_join (
colData (mae) |> as_tibble () |> select (pid, site, arm, ITT, mITT, PP) |> distinct (),
by = join_by (pid, site)
) |>
filter (mITT) |>
select (pid, site, arm, age, race, ethn, cut_size_meals, eat_less, hungry_did_not_eat, sexual_partners_lifetime, sexual_partners_past_month, contraception_method) |>
mutate (
food_insecurity = ifelse (
cut_size_meals == "Yes" | eat_less == "Yes" | hungry_did_not_eat == "Yes" ,
"Yes" ,
"No"
),
food_insecurity = factor (food_insecurity, levels = c ("No" , "Yes" ))
) |>
select (- c (eat_less, cut_size_meals, hungry_did_not_eat))
Code
table1_data |>
ggplot (aes (x = food_insecurity)) +
geom_bar () +
labs (x= "Food insecurity" , y = "Count" )+
theme_classic ()
We also need to include the participants’ partner’s gender.
Participants are asked about their partner’s gender at visits 1000, 1100, 1200, 1300, 1400, 1500, 1700, 1900, 2120 in CRF 19.
The variable past_week_sex_partner of the mae@metadata$visits_crfs_merged table will be summarize into one of these 4 categories:
No sex
Only male
Only female
Both
We add this variable to the table1_data table.
Code
gender_sexual_partner <-
metadata (mae)$ visits_crfs_merged |>
select (pid, past_week_sex_partner) |>
filter (! is.na (past_week_sex_partner)) |>
group_by (pid) |>
summarise (
gender_sexual_partner =
case_when (
all (past_week_sex_partner == "No sex in the past week" ) ~ "No sex" ,
any (past_week_sex_partner %in% c ("A single male partner" , "Multiple male partners" )) &
! any (past_week_sex_partner %in% c ("A single female partner" , "Multiple female partners" )) ~ "Only man" ,
any (past_week_sex_partner %in% c ("A single female partner" , "Multiple female partners" )) &
! any (past_week_sex_partner %in% c ("A single male partner" , "Multiple male partners" )) ~ "Only female" ,
any (past_week_sex_partner %in% c ("A single male partner" , "Multiple male partners" )) &
any (past_week_sex_partner %in% c ("A single female partner" , "Multiple female partners" )) ~ "Both" ,
TRUE ~ NA_character_
),
.groups = "drop"
)
table1_data <-
table1_data |>
left_join (gender_sexual_partner, by = "pid" )
We also re-code the race variable as follows:
“Asian or South Asian” is re-coded to “Asian” (for brievity)
“Black, African American, or African” and “African” are re-coded to “Black” (for brievity)
“Coloured”, “White”, “Other”, and “Prefered not to answer” are left as is.
Code
table1_data <-
table1_data |>
mutate (
across (where (is.character), as.factor),
race = fct_recode (
race,
"Asian" = "Asian or South Asian" ,
"Black" = "Black, African American, or African" ,
"Black" = "African" ,
"Coloured" = "Coloured" ,
"White" = "White" ,
"Other" = "Other" ,
"Prefer not to answer" = "Prefer not to answer"
),
race = fct_relevel (race, "Asian" , "Black" , "Coloured" , "White" , "Other" , "Prefer not to answer" ),
arm = arm |> fct_drop ()
) |>
set_variable_labels (
food_insecurity = "Report food insecurity in past 12 months " ,
sexual_partners_past_month = "Number of partners in past month" ,
sexual_partners_lifetime = "Number of partners in lifetime" ,
ethn = "Ethnicity" ,
contraception_method = "Contraception" ,
gender_sexual_partner = "Gender of sexual partners"
)
Create Table 1 (by arm)
Code
table1_data |>
select (- pid) |>
rename_with (~ str_to_title (.x)) |>
tbl_summary (
by = Arm,
type =
list (
Sexual_partners_past_month ~ "continuous" ,
Sexual_partners_lifetime ~ "continuous"
),
statistic = list (
all_continuous () ~ "{mean} ({min}-{max})" ,
all_categorical () ~ "{n} ({p}%)"
)
) |>
add_p (include = - Site)
Site
CAP
14 (74%)
14 (67%)
14 (93%)
14 (93%)
14 (70%)
MGH
5 (26%)
7 (33%)
1 (6.7%)
1 (6.7%)
6 (30%)
Age
29 (20-38)
29 (18-40)
26 (21-34)
26 (19-35)
29 (20-40)
0.5
Race
0.5
Asian
1 (5.3%)
0 (0%)
1 (6.7%)
0 (0%)
0 (0%)
Black
14 (74%)
18 (86%)
14 (93%)
15 (100%)
18 (90%)
Coloured
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
White
2 (11%)
3 (14%)
0 (0%)
0 (0%)
1 (5.0%)
Other
1 (5.3%)
0 (0%)
0 (0%)
0 (0%)
1 (5.0%)
Prefer not to answer
1 (5.3%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Ethnicity
0.8
Not hispanic/latino/spanish
3 (60%)
6 (86%)
1 (100%)
1 (100%)
4 (67%)
hispanic/latino/spanish
2 (40%)
1 (14%)
0 (0%)
0 (0%)
2 (33%)
Unknown
14
14
14
14
14
Number of partners in lifetime
6 (1-10)
11 (1-90)
4 (1-10)
6 (1-20)
5 (1-10)
0.8
Number of partners in past month
1 (1-2)
1 (0-2)
1 (0-2)
1 (0-3)
1 (0-2)
0.4
Contraception
>0.9
Continuous combined oral contraceptive pills (skip placebo)
3 (16%)
5 (24%)
1 (6.7%)
1 (6.7%)
3 (15%)
Contraceptive vaginal ring
0 (0%)
1 (4.8%)
0 (0%)
0 (0%)
0 (0%)
Cyclic combined oral contraceptive pills
0 (0%)
0 (0%)
0 (0%)
0 (0%)
1 (5.0%)
Depo provera
11 (58%)
10 (48%)
12 (80%)
11 (73%)
11 (55%)
Hormonal implant
1 (5.3%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Hormone containing IUD
2 (11%)
2 (9.5%)
0 (0%)
0 (0%)
2 (10%)
Net-EN
2 (11%)
2 (9.5%)
2 (13%)
3 (20%)
2 (10%)
Other
0 (0%)
1 (4.8%)
0 (0%)
0 (0%)
1 (5.0%)
Report food insecurity in past 12 months
7 (37%)
7 (33%)
6 (40%)
5 (33%)
6 (30%)
>0.9
Gender of sexual partners
0.8
No sex
1 (5.3%)
3 (14%)
1 (6.7%)
1 (6.7%)
3 (15%)
Only man
18 (95%)
18 (86%)
14 (93%)
14 (93%)
17 (85%)
Create Table 1 (by site)
Code
table1_data |>
select (arm, age, race, site, food_insecurity, contraception_method, sexual_partners_past_month, sexual_partners_lifetime) |>
rename_with (~ str_to_title (.x)) |>
tbl_summary (
by = Site,
type =
list (
Sexual_partners_past_month ~ "continuous" ,
Sexual_partners_lifetime ~ "continuous"
),
statistic = list (
all_continuous () ~ "{mean} ({min}-{max})" ,
all_categorical () ~ "{n} ({p}%)"
)
) |>
add_p ()
Arm
0.2
Pl
14 (20%)
5 (25%)
LC-106-7
14 (20%)
7 (35%)
LC-106-3
14 (20%)
1 (5.0%)
LC-106-o
14 (20%)
1 (5.0%)
LC-115
14 (20%)
6 (30%)
Age
27 (18-37)
32 (22-40)
<0.001
Race
<0.001
Asian
0 (0%)
2 (10%)
Black
70 (100%)
9 (45%)
Coloured
0 (0%)
0 (0%)
White
0 (0%)
6 (30%)
Other
0 (0%)
2 (10%)
Prefer not to answer
0 (0%)
1 (5.0%)
Report food insecurity in past 12 months
28 (40%)
3 (15%)
0.038
Contraception
<0.001
Continuous combined oral contraceptive pills (skip placebo)
5 (7.1%)
8 (40%)
Contraceptive vaginal ring
0 (0%)
1 (5.0%)
Cyclic combined oral contraceptive pills
0 (0%)
1 (5.0%)
Depo provera
54 (77%)
1 (5.0%)
Hormonal implant
0 (0%)
1 (5.0%)
Hormone containing IUD
0 (0%)
6 (30%)
Net-EN
11 (16%)
0 (0%)
Other
0 (0%)
2 (10%)
Number of partners in past month
1 (0-3)
1 (0-2)
0.4
Number of partners in lifetime
4 (1-20)
14 (4-90)
<0.001
And finally we create one Table 1 for each site
Code
table1_data |>
filter (site == "CAP" ) |>
select (- c (pid, site)) |>
rename_with (~ str_to_title (.x)) |>
tbl_summary (
by = Arm,
type =
list (
Sexual_partners_past_month ~ "continuous" ,
Sexual_partners_lifetime ~ "continuous"
),
statistic = list (
all_continuous () ~ "{mean} ({min}-{max})" ,
all_categorical () ~ "{n} ({p}%)"
)
) |>
add_p ()|>
modify_caption ("**Site = CAP**" )
Site = CAP
Age
28 (20-36)
26 (18-36)
26 (21-34)
26 (19-35)
27 (20-37)
>0.9
Race
>0.9
Asian
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Black
14 (100%)
14 (100%)
14 (100%)
14 (100%)
14 (100%)
Coloured
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
White
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Other
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Prefer not to answer
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Ethnicity
Not hispanic/latino/spanish
0 (NA%)
0 (NA%)
0 (NA%)
0 (NA%)
0 (NA%)
hispanic/latino/spanish
0 (NA%)
0 (NA%)
0 (NA%)
0 (NA%)
0 (NA%)
Unknown
14
14
14
14
14
Number of partners in lifetime
5 (1-10)
4 (1-7)
4 (1-10)
6 (1-20)
4 (1-9)
0.8
Number of partners in past month
1 (1-2)
1 (1-2)
1 (0-2)
1 (0-3)
1 (1-2)
0.2
Contraception
0.9
Continuous combined oral contraceptive pills (skip placebo)
1 (7.1%)
2 (14%)
0 (0%)
0 (0%)
2 (14%)
Contraceptive vaginal ring
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Cyclic combined oral contraceptive pills
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Depo provera
11 (79%)
10 (71%)
12 (86%)
11 (79%)
10 (71%)
Hormonal implant
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Hormone containing IUD
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Net-EN
2 (14%)
2 (14%)
2 (14%)
3 (21%)
2 (14%)
Other
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Report food insecurity in past 12 months
7 (50%)
5 (36%)
6 (43%)
5 (36%)
5 (36%)
>0.9
Gender of sexual partners
>0.9
No sex
1 (7.1%)
2 (14%)
1 (7.1%)
1 (7.1%)
2 (14%)
Only man
13 (93%)
12 (86%)
13 (93%)
13 (93%)
12 (86%)
Code
table1_data |>
filter (site == "MGH" ) |>
select (- c (pid, site)) |>
rename_with (~ str_to_title (.x)) |>
tbl_summary (
by = Arm,
type =
list (
Sexual_partners_past_month ~ "continuous" ,
Sexual_partners_lifetime ~ "continuous"
),
statistic = list (
all_continuous () ~ "{mean} ({min}-{max})" ,
all_categorical () ~ "{n} ({p}%)"
)
) |>
add_p ()|>
modify_caption ("**Site = MGH**" )
Site = MGH
Age
32 (26-38)
33 (24-40)
25 (25-25)
35 (35-35)
33 (22-40)
0.7
Race
0.13
Asian
1 (20%)
0 (0%)
1 (100%)
0 (0%)
0 (0%)
Black
0 (0%)
4 (57%)
0 (0%)
1 (100%)
4 (67%)
Coloured
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
White
2 (40%)
3 (43%)
0 (0%)
0 (0%)
1 (17%)
Other
1 (20%)
0 (0%)
0 (0%)
0 (0%)
1 (17%)
Prefer not to answer
1 (20%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Ethnicity
0.8
Not hispanic/latino/spanish
3 (60%)
6 (86%)
1 (100%)
1 (100%)
4 (67%)
hispanic/latino/spanish
2 (40%)
1 (14%)
0 (0%)
0 (0%)
2 (33%)
Number of partners in lifetime
9 (7-10)
25 (6-90)
7 (7-7)
10 (10-10)
7 (4-10)
0.059
Number of partners in past month
1 (1-2)
1 (0-2)
1 (1-1)
1 (1-1)
1 (0-2)
>0.9
Contraception
>0.9
Continuous combined oral contraceptive pills (skip placebo)
2 (40%)
3 (43%)
1 (100%)
1 (100%)
1 (17%)
Contraceptive vaginal ring
0 (0%)
1 (14%)
0 (0%)
0 (0%)
0 (0%)
Cyclic combined oral contraceptive pills
0 (0%)
0 (0%)
0 (0%)
0 (0%)
1 (17%)
Depo provera
0 (0%)
0 (0%)
0 (0%)
0 (0%)
1 (17%)
Hormonal implant
1 (20%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Hormone containing IUD
2 (40%)
2 (29%)
0 (0%)
0 (0%)
2 (33%)
Net-EN
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
Other
0 (0%)
1 (14%)
0 (0%)
0 (0%)
1 (17%)
Report food insecurity in past 12 months
0 (0%)
2 (29%)
0 (0%)
0 (0%)
1 (17%)
0.8
Gender of sexual partners
>0.9
No sex
0 (0%)
1 (14%)
0 (0%)
0 (0%)
1 (17%)
Only man
5 (100%)
6 (86%)
1 (100%)
1 (100%)
5 (83%)
Primary outcomes (Table 2)
Colonization by week 5
Code
col_by_week5 <-
mae@ colData |> as_tibble () |>
select (pid, site, randomized, arm, ITT, mITT, PP) |>
distinct () |>
filter (randomized, mITT) |>
left_join (
mae[["col_LBP_mg" ]] |>
as_tibble () |>
dplyr:: left_join (
mae@ colData |> as_tibble () |> select (uid, pid, visit_code),
by = join_by (.sample == uid)
) |>
dplyr:: filter (visit_code == "1500" , .feature == "colonized_LBP_by_mg" )
) |>
mutate (LBP_colonization_by_week5 = outcome |> replace_na (FALSE ))
# by qPCR
col_by_week5_qPCR <-
mae@ colData |> as_tibble () |>
select (pid, site, randomized, arm, ITT, mITT, PP) |>
distinct () |>
filter (randomized, mITT) |>
left_join (
mae[["col_LBP_qPCR" ]] |>
as_tibble () |>
dplyr:: left_join (
mae@ colData |> as_tibble () |> select (uid, pid, visit_code),
by = join_by (.sample == uid)
) |>
dplyr:: filter (visit_code == "1500" , .feature == "colonized_LBP_by_qpcr" )
) |>
mutate (LBP_colonization_by_week5 = outcome |> replace_na (FALSE ))
detected_by_week5_qPCR <-
mae@ colData |> as_tibble () |>
select (pid, site, randomized, arm, ITT, mITT, PP) |>
distinct () |>
filter (randomized, mITT) |>
left_join (
mae[["col_LBP_qPCR" ]] |>
as_tibble () |>
dplyr:: left_join (
mae@ colData |> as_tibble () |> select (uid, pid, visit_code),
by = join_by (.sample == uid)
) |>
dplyr:: filter (visit_code == "1500" , .feature == "LBP_detected_by_qpcr" )
) |>
mutate (LBP_detected_by_week5 = outcome |> replace_na (FALSE ))
# crispatus dominance
crisp_dom_by_5week <-
mae@ colData |> as_tibble () |>
select (pid, site, randomized, arm, ITT, mITT, PP) |>
distinct () |>
filter (randomized, mITT) |>
left_join (
mae[["col_crispatus_mg" ]] |>
as_tibble () |>
dplyr:: left_join (
mae@ colData |> as_tibble () |> select (uid, pid, visit_code),
by = join_by (.sample == uid)
) |>
dplyr:: filter (visit_code == "1500" , .feature == "crispatus_dominance_by_mg" )
) |>
mutate (crisp_dom_by_5week = outcome |> replace_na (FALSE ))
Visualization of the primary outcome (and a few selected secondary outcomes)
Code
col_by_week5 |>
arrange (site, arm, LBP_colonization_by_week5) |>
group_by (site, arm) |>
mutate (participant_number = row_number () |> factor ()) |>
ungroup () |>
ggplot (aes (x = participant_number, y = site |> fct_rev (), col = site, fill = site, shape = LBP_colonization_by_week5)) +
geom_point (size = 3 ) +
facet_grid (. ~ arm, scales = "free_x" , space = "free_x" ) +
guides (color = "none" , fill = "none" ) +
scale_shape_manual ("Colonized by LBP by week 5" , values = c (1 , 16 )) +
scale_color_manual (values = site_colors) +
xlab ("Participant number" ) + ylab ("" ) +
theme (
legend.position = "top" ,
axis.text.x = element_text (angle = 90 , hjust = 1 , vjust = 0.5 )
)
Code
col_by_week5_qPCR |>
arrange (site, arm, LBP_colonization_by_week5) |>
group_by (site, arm) |>
mutate (participant_number = row_number () |> factor ()) |>
ungroup () |>
ggplot (aes (x = participant_number, y = site |> fct_rev (), col = site, fill = site, shape = LBP_colonization_by_week5)) +
geom_point (size = 3 ) +
facet_grid (. ~ arm, scales = "free_x" , space = "free_x" ) +
guides (color = "none" , fill = "none" ) +
scale_shape_manual ("Colonized by LBP by week 5" , values = c (1 , 16 )) +
scale_color_manual (values = site_colors) +
xlab ("Participant number" ) + ylab ("" ) +
ggtitle ("Colonization based on the relative abundances estimated by qPCR" ) +
theme (
legend.position = "top" ,
axis.text.x = element_text (angle = 90 , hjust = 1 , vjust = 0.5 ),
plot.title = element_text (hjust = 0.5 )
)
Code
detected_by_week5_qPCR |>
arrange (site, arm, LBP_detected_by_week5) |>
group_by (site, arm) |>
mutate (participant_number = row_number () |> factor ()) |>
ungroup () |>
ggplot (aes (x = participant_number, y = site |> fct_rev (), col = site, fill = site, shape = LBP_detected_by_week5)) +
geom_point (size = 3 ) +
facet_grid (. ~ arm, scales = "free_x" , space = "free_x" ) +
guides (color = "none" , fill = "none" ) +
scale_shape_manual ("Colonized by LBP by week 5" , values = c (1 , 16 )) +
scale_color_manual (values = site_colors) +
xlab ("Participant number" ) + ylab ("" ) +
ggtitle ("At least 2 LBP strains detected at the same visit by week 5 qPCR" ) +
theme (
legend.position = "top" ,
axis.text.x = element_text (angle = 90 , hjust = 1 , vjust = 0.5 ),
plot.title = element_text (hjust = 0.5 )
)
Code
crisp_dom_by_5week |>
arrange (site, arm, crisp_dom_by_5week) |>
group_by (site, arm) |>
mutate (participant_number = row_number () |> factor ()) |>
ungroup () |>
ggplot (aes (x = participant_number, y = site |> fct_rev (), col = site, fill = site, shape = crisp_dom_by_5week)) +
geom_point (size = 3 ) +
facet_grid (. ~ arm, scales = "free_x" , space = "free_x" ) +
guides (color = "none" , fill = "none" ) +
scale_shape_manual ("Colonized by LBP by week 5" , values = c (1 , 16 )) +
scale_color_manual (values = site_colors) +
xlab ("Participant number" ) + ylab ("" ) +
ggtitle ("L. crispatus dominance (≥ 50% by metagenomics)" ) +
theme (
legend.position = "top" ,
axis.text.x = element_text (angle = 90 , hjust = 1 , vjust = 0.5 ),
plot.title = element_text (hjust = 0.5 )
)
Code
summary_table <-
col_by_week5 |>
group_by (site, arm) |>
summarize (
n = n (),
n_success = sum (LBP_colonization_by_week5),
.groups = "drop"
) |>
mutate (
p = n_success / n,
CI = binom:: binom.confint (n_success, n, method = "exact" )
)
Code
table_2 <-
summary_table |>
mutate (
N_per_site_and_arm = str_c ("N = " , n),
LBP_strain_detected = str_c (n_success, " (" , round (p * 100 )," %)" ),
CI_text = str_c (round (CI$ lower * 100 ), "% - " , round (CI$ upper * 100 ), "%" )
) |>
dplyr:: select (
site, arm, N_per_site_and_arm, LBP_strain_detected, CI_text
) |>
pivot_longer (- c (site, arm)) |>
mutate (
name =
case_when (
name == "N_per_site_and_arm" ~ "N participants" ,
name == "LBP_strain_detected" ~ "n (%) participants \n with LBP strain detected" ,
name == "CI_text" ~ "95% CI"
)
) |>
pivot_wider (
id_cols = c (site, name),
names_from = arm,
values_from = value, values_fill = ""
)
Code
table_2 |>
group_by (site) |>
gt (caption = "Table 1: Colonization by week 5 by arm and site" , row_group_as_column = TRUE ) |>
cols_width (
"name" ~ px (200 ),
everything () ~ px (120 )
) |>
cols_label (
name = "" ,
# Blinded = "All blinded arms",
Pl = "Placebo" ,
` LC-106-7 ` = "LC-106<br>7 days" ,
` LC-106-3 ` = "LC-106<br>3 days" ,
` LC-106-o ` = "LC-106<br>overlap" ,
` LC-115 ` = "LC-115<br>7 days" ,
.fn = md
)
Table 1: Colonization by week 5 by arm and site
CAP
N participants
N = 14
N = 14
N = 14
N = 14
N = 14
n (%) participants with LBP strain detected
0 (0 %)
9 (64 %)
7 (50 %)
7 (50 %)
10 (71 %)
95% CI
0% - 23%
35% - 87%
23% - 77%
23% - 77%
42% - 92%
MGH
N participants
N = 5
N = 7
N = 1
N = 1
N = 6
n (%) participants with LBP strain detected
1 (20 %)
6 (86 %)
1 (100 %)
1 (100 %)
6 (100 %)
95% CI
1% - 72%
42% - 100%
3% - 100%
3% - 100%
54% - 100%
Tests
We test for differences between the placebo arm and each active arm using Fisher’s exact test (exhaustive permutation test).
Data from both site are merged for each arms given the low numbers of participants at MGH.
The \(p\) -values are adjusted for multiple testing using the Benjamini-Hochberg method.
Code
arms_to_test <- col_by_week5$ arm |> unique () |> setdiff ("Pl" )
fisher_test_results <-
map (
arms_to_test,
~ {
arm_to_test <- .x
tmp <- col_by_week5 |> filter (arm %in% c ("Pl" , arm_to_test))
table_for_test <- table (tmp$ LBP_colonization_by_week5, tmp$ arm |> fct_drop ())
test_res <- fisher.test (table_for_test, alternative = "greater" )
tibble (
arm = arm_to_test,
p_value = test_res$ p.value,
OR = test_res$ estimate,
lower_CI = test_res$ conf.int[1 ],
upper_CI = test_res$ conf.int[2 ]
)
}
) |>
bind_rows () |>
mutate (
adj_p_value = p_value |> p.adjust (method = "BH" )
)
Code
test_results <-
fisher_test_results |>
mutate (
adj_p_value = adj_p_value |> scales:: pvalue (accuracy = 0.001 ),
OR = OR |> round (2 ) |> as.character (),
CI = str_c (lower_CI |> round (2 ), " - " , upper_CI |> round (2 ))
) |>
select (arm, adj_p_value, OR, CI) |>
pivot_longer (cols = - arm, names_to = "name" , values_to = "value" ) |>
pivot_wider (names_from = arm, values_from = value) |>
mutate (
name =
case_when (
name == "adj_p_value" ~ "Adj. p-value" ,
name == "OR" ~ "OR" ,
name == "CI" ~ "95% CI" ,
TRUE ~ name
)
)
Code
table_2_with_model_results <-
bind_rows (
table_2,
test_results
) |>
mutate (
site = site |> str_replace_na ("" ),
Pl =
case_when (
name == "Adj. p-value" ~ "Ref." ,
name == "OR" ~ "" ,
str_detect (name, "CI" ) ~ "" ,
TRUE ~ Pl
)
)
Code
table_2_with_model_results |>
group_by (site) |>
gt (caption = "Table 2: Colonization by week 5 by arm and site" , row_group_as_column = TRUE ) |>
cols_width (
"name" ~ px (150 ),
everything () ~ px (120 )
) |>
cols_label (
name = "" ,
Pl = "Placebo" ,
` LC-106-7 ` = "LC-106<br>7 days" ,
` LC-106-3 ` = "LC-106<br>3 days" ,
` LC-106-o ` = "LC-106<br>overlap" ,
` LC-115 ` = "LC-115<br>7 days" ,
.fn = md
)
Table 2: Colonization by week 5 by arm and site
CAP
N participants
N = 14
N = 14
N = 14
N = 14
N = 14
n (%) participants with LBP strain detected
0 (0 %)
9 (64 %)
7 (50 %)
7 (50 %)
10 (71 %)
95% CI
35% - 87%
23% - 77%
23% - 77%
42% - 92%
MGH
N participants
N = 5
N = 7
N = 1
N = 1
N = 6
n (%) participants with LBP strain detected
1 (20 %)
6 (86 %)
1 (100 %)
1 (100 %)
6 (100 %)
95% CI
42% - 100%
3% - 100%
3% - 100%
54% - 100%
Adj. p-value
Ref.
<0.001
0.002
0.002
<0.001
OR
39.76
18.61
18.61
60.46
95% CI
5.69 - Inf
2.48 - Inf
2.48 - Inf
8.17 - Inf
Model
Code
col_by_week5_agg <-
col_by_week5 |>
group_by (site, arm) |>
summarise (success = LBP_colonization_by_week5 |> sum (), ntrials = n ())
bfit <-
brm (
success | trials (ntrials) ~ arm + site + site: arm,
data = col_by_week5_agg,
family = binomial ("logit" )
)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 2.6e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.81 seconds (Warm-up)
Chain 1: 0.884 seconds (Sampling)
Chain 1: 1.694 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 4e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.734 seconds (Warm-up)
Chain 2: 0.973 seconds (Sampling)
Chain 2: 1.707 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 4e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.834 seconds (Warm-up)
Chain 3: 0.864 seconds (Sampling)
Chain 3: 1.698 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.935 seconds (Warm-up)
Chain 4: 0.951 seconds (Sampling)
Chain 4: 1.886 seconds (Total)
Chain 4:
Code
summary (bfit)
plot (bfit, ask = FALSE )
bind_cols (col_by_week5_agg, predict (bfit))
# loo(bfit, moment_match = TRUE, reloo = TRUE)
# pp_check(bfit)
Code
plot_res <- plot (conditional_effects (bfit), ask = FALSE )
Code
plot_res <- lapply (plot_res, function (p) {
p + scale_color_manual (values = site_colors)
})
Reduce (` + ` , plot_res) +
plot_layout (widths = c (2.2 , 1 , 3 ))
Code
model_results <- summary (bfit)$ fixed
model_results <-
model_results |>
rownames_to_column ("term" ) |>
as_tibble () |>
mutate (
factor =
case_when (
str_detect (term, ":site" ) ~ "Arm x Site" ,
str_detect (term, "^site" ) ~ "Site" ,
str_detect (term, "^arm" ) ~ "Arm" ,
TRUE ~ "Ref." ,
),
factor_level =
term |>
str_remove ("arm" ) |> str_remove ("site" ) |> str_replace (":" ," x " ) |>
str_replace ("6M" , "6-" ) |> str_replace ("CM" , "C-" )
) |>
mutate (
OR = ifelse (factor == "Ref." , NA , exp (Estimate)),
lower_CI = ifelse (factor == "Ref." , NA , exp (` l-95% CI ` )),
upper_CI = ifelse (factor == "Ref." , NA , exp (` u-95% CI ` ))
)
Code
model_results |>
select (factor, factor_level, OR, lower_CI, upper_CI) |>
gt () |>
gt:: fmt_number (
decimals = 2 , n_sigfig = 3
)
Code
model_results_arm <-
model_results |>
filter (factor == "Arm" ) |>
mutate (
OR = OR |> round (2 ) |> as.character (),
CI = str_c (lower_CI |> round (2 ), " - " , upper_CI |> round (2 ))
) |>
select (factor_level, OR, CI) |>
pivot_longer (cols = - factor_level, names_to = "name" , values_to = "value" ) |>
pivot_wider (names_from = factor_level, values_from = value)
Code
table_2_with_model_results <-
bind_rows (
table_2,
model_results_arm
) |>
mutate (
site = site |> str_replace_na ("" ),
Pl =
case_when (
name == "OR" ~ "Ref." ,
name == "CI" ~ "" ,
TRUE ~ Pl
)
)
Code
table_2_with_model_results |>
group_by (site) |>
gt (caption = "Table 2: Colonization by week 5 by arm and site" , row_group_as_column = TRUE ) |>
cols_width (
"name" ~ px (200 ),
everything () ~ px (120 )
) |>
cols_label (
name = "" ,
Pl = "Placebo" ,
` LC-106-7 ` = "LC-106<br>7 days" ,
` LC-106-3 ` = "LC-106<br>3 days" ,
` LC-106-o ` = "LC-106<br>overlap" ,
` LC-115 ` = "LC-115<br>7 days" ,
.fn = md
)